Week 2: Data Visualization

library(tidyverse) #ggplot, readr, dplyer, etc.
library(haven) #opening Stata datasets
library(skimr) #descriptive statistics
library(nycflights13) #load datasets from moderndive book

states <- read_dta("states.dta")

#install.packages("scales")
library(scales) #for labeling axes with %, $, etc.
#install.packages("ggridges")
library(ggridges)
#install.packages("moderndive")
library(moderndive)
#install.packages("ggrepel")
library(ggrepel)
#install.packages("questionr")
library(questionr) #for frequency tables

nes <- read_dta("nes.dta")
states <- read_dta("states.dta")
vdem <- read_dta("vdem.dta")

# Calling up variables in datasets: $ method, "data$variable"
skim(vdem$v2x_polyarchy)

# For skimr, you can also use use (data, variable)
skim(vdem, v2x_polyarchy)

# Frequency table, using questionr package: 
freq(vdem$v2x_regime)

# We'll use the "weather" data from nycflights13 package; 
# Hourly meteorological data for LGA, JFK and EWR.
# We want to visualize the distribution of a single variable with a histogram; 
# these are ideal for "continuous" variables with lots of values.
weather <- weather

# We'll focus on temperatures ("temp") at these three airports
# Note: geometry for a histogram is "geom_histogram"
ggplot(data = weather, mapping = aes(x = temp)) +
  geom_histogram() 

# Too clumpy; let's delineate the bars some more. Always use color = "white" option.
ggplot(data = weather, mapping = aes(x = temp)) +
  geom_histogram(color = "white")

# We can also adjust the "bins" or the width of the bars; more or less fine-grained
ggplot(data = weather, mapping = aes(x = temp)) +
  geom_histogram(bins = 40, color = "white")

ggplot(data = weather, mapping = aes(x = temp)) +
  geom_histogram(binwidth = 10, color = "white")

# Rules for binwidth? The Freedman-Diaconis rule:
# binwwidth <- 2 * IQR / (N^(1/3))
skim(weather, temp)
bw <- (2 * (70.0-39.9)) / (26115^(1/3))

ggplot(data = weather, mapping = aes(x = temp)) +
  geom_histogram(binwidth = bw, color = "white")

# Facet wraps: These are very cool for subsetting by some grouping

ggplot(data = weather, mapping = aes(x = temp)) +
  geom_histogram(binwidth = 5, color = "white") +
  facet_wrap(~month)

# We can also report y-axis in percentage terms; the third line uses 
# the "scales" package
# also, label axes, give title
ggplot(data = weather, mapping = aes(x = temp)) +
  geom_histogram(aes(y = (..count..)/sum(..count..)), binwidth = 2.029083, color = "white") +
  scale_y_continuous(labels=percent) +
  labs(x="Temperature", y="Percentage of Hours", title="Temperatures")

# Now change labeling for y-axis values
ggplot(data = weather, mapping = aes(x = temp)) +
  geom_histogram(aes(y = (..count..)/sum(..count..)), binwidth = 2.029083, color = "white") +
  scale_y_continuous(limits = c(0, .06), breaks = seq(0, .06, by = .01), labels=percent) +
  labs(x="Temperature", y="Percentage of Hours", title="Temperatures")

# Change theme for graph background
# Now change labeling for y-axis values
ggplot(data = weather, mapping = aes(x = temp)) +
  geom_histogram(aes(y = (..count..)/sum(..count..)), binwidth = 2.029083, color = "white") +
  scale_y_continuous(limits = c(0, .06), breaks = seq(0, .06, by = .01), labels=percent) +
  labs(x="Temperature", y="Percentage of Hours", title="Temperatures") +
  theme_minimal()

# Center title
ggplot(data = weather, mapping = aes(x = temp)) +
  geom_histogram(aes(y = (..count..)/sum(..count..)), binwidth = 2.029083, color = "white") +
  scale_y_continuous(limits = c(0, .06), breaks = seq(0, .06, by = .01), labels=percent) +
  labs(x="Temperature", y="Percentage of Hours", title="Temperatures") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

# There are some aspects of the syntax used above that are unnecessary and that we can
# eliminate. For example: 

# This command: 
ggplot(data = weather, mapping = aes(x = temp)) +
  geom_histogram()

# is identical to this commmand (note that I've removed "data=" and "mapping="): 
ggplot(weather, aes(x = temp)) +
  geom_histogram()

# You actually don't even need "x=" but I like to include that to make it clear.
ggplot(weather, aes(temp)) +
  geom_histogram()

# From now on, we'll simplify to the second command

ggplot(weather, aes(x = temp)) +
  geom_density() +
  labs(x="Temperature", y="Density", title="Temperatures") +
  theme(plot.title = element_text(hjust = 0.5))

# Fill with color
ggplot(weather, aes(x = temp)) +
  geom_density(fill="dodgerblue") +
  labs(x="Temperature", y="Density", title="Temperatures") +
  theme(plot.title = element_text(hjust = 0.5))

# Make transparent with "alpha"
ggplot(weather, aes(x = temp)) +
  geom_density(fill="dodgerblue", alpha=.5) +
  labs(x="Temperature", y="Density", title="Temperatures") +
  theme(plot.title = element_text(hjust = 0.5))

# Ridge plots
# Notice I have to change month to "factor" variable. Most variables are "continuous"
# by default. I need to change this to "discrete" or ggplot can't graph it like 
# I want it to.
ggplot(weather, aes(x = temp, y=factor(month))) +
  geom_density_ridges(fill="dodgerblue", alpha=.5) +
  labs(x="Temperature", y="Month", title="Temperatures") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

# What do boxplots give us exactly? "Low", 25th pctile, 50th pctile (median), 75th pctile, 
# and "high"; and outliers.

# Let's again use the weather data: temperature 
# Geometry for boxplot is "geom_boxplot"
ggplot(weather, aes(y = temp)) +
  geom_boxplot() 

# Use more fine-grained y-axis labels
ggplot(weather, aes(y = temp)) +
  geom_boxplot() + 
  scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, by = 5))

# Now let's look at temp by month.
# We need to treat month as a "factor" or categorical variable instead of a continuous 
# variable. We'll talk more about this distinction later in the class.
ggplot(weather, aes(x = factor(month), y = temp)) +
  geom_boxplot() + 
  labs(x="Month", y="Temperature")

# Call up flights data
flights <- flights

# We could get a frequency distribution in table form using the 
# "freq" command from the questionr package.
freq(flights$carrier)

# Now ggplot; geometry for bar graph is "geom_bar"
ggplot(flights, aes(x = carrier)) +
  geom_bar()

# Stacked bar graphs; let's say you wanted to break down by origin (EWR, JFK, or LGA)
ggplot(flights, aes(x = carrier, fill = origin)) +
  geom_bar()

# Grouped bar graphs; let's say you wanted to break down by origin (EWR, JFK, or LGA)
ggplot(flights, aes(x = carrier, fill = origin)) +
  geom_bar(position = "dodge")

# We could also facet wrap
ggplot(flights, aes(x = carrier)) +
  geom_bar() +
  facet_wrap(~ origin, ncol = 1)

# Note how y-axis is a count by default. We could report as a proportion instead.
ggplot(flights, aes(x = carrier)) +
  geom_bar(aes(y = (..count..)/sum(..count..)))

# Percentage much more intuitive; the third line uses the package "scales", which 
# converts the proportion to a percentage and adds a percent sign.
ggplot(flights, aes(x = carrier)) +
  geom_bar(aes(y = (..count..)/sum(..count..))) +
  scale_y_continuous(labels=percent) 

# Report y-axis labels as percentage; the third line uses the package "scales"; label axes
ggplot(flights, aes(x = carrier)) +
  geom_bar(aes(y = (..count..)/sum(..count..))) +
  scale_y_continuous(labels=percent) +
  labs(x="Carrier", y="Percentage")

# Alaska data (from moderndive package); let's view it first
alaska_flights <- alaska_flights

# Variable descriptions
#?alaska_flights

# Basic scatterplot; note how the plus sign separates different arguments.
# Note how the three elements of "grammar of graphics" are implemented with ggplot.
# The geometry name for a scatterplot is "geom_point"
ggplot(alaska_flights, aes(x = dep_delay, y = arr_delay)) + 
  geom_point()

# We can add different features, which we'll do throughout the semester. 
# Let's add more informative labels for y and x axis
ggplot(alaska_flights, aes(x = dep_delay, y = arr_delay)) + 
  geom_point() +
  labs(x="Departure Delay", y="Arrival Delay")

# Transparent dots and white background
ggplot(alaska_flights, aes(x = dep_delay, y = arr_delay)) + 
  geom_point(alpha = .2) +
  labs(x="Departure Delay", y="Arrival Delay") +
  theme_minimal()

# Make dots smaller
ggplot(alaska_flights, aes(x = dep_delay, y = arr_delay)) + 
  geom_point(size=.5, alpha=.2) +
  labs(x="Departure Delay", y="Arrival Delay") +
  theme_minimal()

# "Jitter" the points -- good option if there are many overlaps
ggplot(alaska_flights, aes(x = dep_delay, y = arr_delay)) + 
  geom_jitter(width = 30, height = 30) +
  labs(x="Departure Delay", y="Arrival Delay") +
  theme_minimal()

# Color code by month
ggplot(alaska_flights, aes(x = dep_delay, y = arr_delay, color=factor(month))) +
  geom_jitter(width = 30, height = 30) +
  labs(x="Departure Delay", y="Arrival Delay") +
  theme_minimal()

# Let's add a "line of best fit" to the plot
ggplot(alaska_flights, aes(x = dep_delay, y = arr_delay)) + 
  geom_point() +
  geom_smooth(method="lm") +
  labs(x="Departure Delay", y="Arrival Delay") +
  theme_minimal()

# Take out the confidence interval (shading around line)
ggplot(alaska_flights, aes(x = dep_delay, y = arr_delay)) + 
  geom_point() +
  geom_smooth(method="lm", se=FALSE) +
  labs(x="Departure Delay", y="Arrival Delay") +
  theme_minimal()

# Set limits and increments for y-axis labels
ggplot(alaska_flights, aes(x = dep_delay, y = arr_delay)) + 
  geom_point() +
  geom_smooth(method="lm", se=FALSE) +
  labs(x="Departure Delay", y="Arrival Delay") +
  scale_y_continuous(limits = c(-100, 200), breaks = seq(-100, 200, by = 50)) +
  theme_minimal()

# We'll use the "early_january_weather" data from moderndive package
early_january_weather <- early_january_weather
#?early_january_weather

# Let's plot temperatures ("temp") in Jan. against time ("time_hour")
# The geometry for a linegraph is "geom_line"
ggplot(early_january_weather, aes(x = time_hour, y = temp)) +
  geom_line()

# Smoother to envision trend and other options. 
ggplot(early_january_weather, aes(x = time_hour, y = temp)) +
  geom_line() +
  geom_smooth(se=FALSE) +
  labs(x="Time/hour", y="Temperature") +
  scale_y_continuous(limits = c(0, 75), breaks = seq(0, 75, by = 10))